Mississippi State University 
Department of Aerospace Engineering 
Drawer A 

Mississippi State, MS 39762 
601-325-3623 


Improved 3~D Turbomachinery CFD Algorithm 
November 1, 1986 - May 15, 1988 


June 1988 


Final Report 

J. Mark Janus 
David L. tyhitfield 


Prepared for 


Scientific Research Associates, Inc. 
Glastonbury, CT 06033 


(NASA-CR-1 793 64) IMPROVED 3-D 
TUB BOM ACHINERY CFD ALGORITHM Final Report, 1 
NOV./1986 - 15 May 1988 (Mississippi State 
Uni v. ) 56 p CSCL 20D 

G3/34 


N 88- 266 19 

Unclas 

0154098 


ABSTRACT 


This report describes the building blocks of a computer algorithm 
developed for the time-accurate flow analysis of rotating machines. The 
flow model is a finite volume method utilizing a high resolution approx- 

t 

imate Riemann solver for interface flux definitions. This block LU 
implicit numerical scheme possesses apparent unconditional stability. 
Multi-block composite gridding is used to orderly partition the field 
into a specified arrangement. Block interfaces, including dynamic in- 
terfaces, are treated such as to mimic interior block communication. 
Special attention is given to the reduction of in-core memory require- 
ments by placing the burden on secondary storage media. Broad applica- 
bility is implied, although the results presented are restricted to that 
of an even blade count configuration. Several other configurations are 
presently under investigation, the results of which will appear in sub- 


sequent publications. 



INTRODUCTION 


Current interest in the design of highly efficient turbomachinery 
has spurred the development of computer algorithms capable of accurately 
analyzing the flowfields of rotating geometries. These new algorithms 
can be used to complement the present-day algorithms which are based 
primarily on empiricism. In addition, the design process can be accel- 
erated with the help of computational fluid dynamics (CFD) as a prelimi- 
nary geometry analyzer. CFD is progressing to the stage where routine 
initial geometry flow analysis of relatively complex rotating machines 
will be performed to reduce the overall cost of project analysis. This 
is being accomplished by the recent, significant improvements in both 
computer hardware and software. This report is concerned with the soft- 
ware aspect of computational analysis. 

At present, there exists a need for a robust computational algo- 
rithm to accurately simulate the fluid flow within a three-dimensional 
rotating geometrical environment. Of key interest will be computational 
efficiency, solution accuracy, and logical simplicity. Previously, the 
focus of this project was on the computational efficiency of the algo- 
rithm . 1 This was primarily due to a realization that the flowfield in 
question would be so computationally exhaustive that even state-of-the- 
art computational resources would require excessive turn-around time. 
Although algorithm efficiency is still a premier issue, it was evident 
at the onset of the continuation of this project that the present level 
of efficiency must be accepted for the time being. The time had come to 
take the numerical algorithm and incorporate it into a code whose frame- 


work centered around turbomachinery . 



This report outlines the mathematical background material and the 
fundamental code logic necessary for the development of a three-dimen- 
sional computer code suitable for the time-accurate analysis of turbo- 

2 

machinery. A detailed report on this subject matter is currently in 
preparation and not available at the time of this printing. 

In the next section a brief description is given of the mathemati- 
cal equations (tools) used in this effort to develop a flowfield model. 
The third section outlines the procedure to build a coherent mathemati- 
cal model within which each equation is exploited for its beneficial 
attributes. The third section also contains a description of the numer- 
ical scheme developed to solve these modeling equations. The fourth 
section presents the logical approach developed for simulating a dynamic 
cylindrical environment in addition to outlining specific routines de- 
veloped to optimize (tailor) the computer code for use on a vector proc- 
essor, namely a Cray X-MP. The results of a complex rotating 
configuration are presented in section five, demonstrating the utility 
of the turbomachinery computer algorithm. The final section is reserved 
for the conclusions of this entire endeavor along with future plans to 
enhance the generality and solution accuracy of the code. 
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MODELING MATHEMATICS 


Conservative Integral Model 


The most fundamental form of mathematical model describing the be- 
havior of a continuum fluid can be derived using a control volume 
approach in an Eulerian reference frame. For computational brevity an 
assumption of a nonconducting, inviscid, perfect gas with no body forces 
will be made here, yielding 
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where V is a control volume bounded by surface S. 

Using divergence theory and additional vector notation, Eqs. (1) 
can be cast in the following condensed form. 
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where in a Cartesian reference frame 
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For use within the context of this report (non-constant V), the time 
differentiation appearing in Eq. (2) must be performed with the aid of 
Leibnitz's Rule. With the aforementioned assumptions, Eq. (2) is com- 
monly referred to as the Euler equations. 


Conservative Differential Form 

Equation (2) holds for an arbitrary fixed control volume V, conse- 
quently the integrand satisfies 
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This equation, referred to as the differential conservation law form of 
the Euler equations, will prove to be very beneficial in matching the 
numerics with the physics being modeled. 

For generality it is desirable to transform the equations from a 
Cartesian reference- frame to a body fitted curvilinear reference frame. 
The general curvilinear axes are defined as follows 

5 = 5(x,y,z,t) , 

n = n(x,y,z,t) , (*0 

4 = C (x , y , z , t ) , 

T = t . 
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Note, the obvious implication that a body fitted reference frame for a 
dynamic geometry is inherently time dependent. 

p O ll 

The details of the transformation are readily available, hence 

they are not repeated here. Transforming the Cartesian equations to the 
more general set of curvilinear equations (dropping vector symbols) 
yields • 
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with the contravar iant velocities 
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the Jacobian of the inverse transformation, i.e. 3(x,y,z)/3(5,n,c) 
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and the metric quantities 
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Quasilinear Form 
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A mathematical analysis of the differential form can yield signifi- 
cant insight into the physics of the flowfield in question. Guided by 
this insight, specific numerical techniques can be employed to more ac- 
curately model the flowfield. This analysis begins with the derivation 
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of yet another form of the equations. The final form necessary in the 
development of this mathematical model is derived from the differential 
form. 

Consider differentiating the vector-valued flux functions F, G, and 
H (homogeneous functions of degree one in Q) appearing in Eq. (5) with 
respect to the dependent variable vector Q. The resulting flux Jaco- 
bians A, B, and C, are written 
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Substituting these expressions in Eq. (5) yields the quasilinear form of 
the Euler equations. It should be noted that the quasilinear form, 
given below, is not in conservation law form. 
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Riemann Problem 


Essential to the investigation of compressible transonic fluids is 
the study of discontinuous solutions. To complete the development of 
the mathematical model, it is necessary to introduce the concept of a 
particular type of initial-value problem for a hyperbolic system of con- 
servation laws. To begin the discussion, a digression to one-dimension 
is in order. Consider for example, the one-dimensional hyperbolic sys- 
tem 


3u 3f 
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where the initial data are 

I u^ for -» < x £ 0 

, 

Up for 0 S x < » 

and u(x,t) is an m component column vector and f(u) is an m component 
vector-valued flux function. The preceding example of the simplest kind 
of discontinuous data (two semi-infinite regions of constant value, sepa- 
rated by a point of discontinuity) poses what is referred to as a 
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Riemann problem. Riemann problems present themselves in varying degrees 
of difficulty, depending on their linearity or lack thereof and on the 
degrees-of-freedom (m) they represent. For a detailed study of the 
Riemann problem and the mathematical theory of shocks (discontinuities) 
the reader is referred to [5] . 
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NUMERICAL FORMULAS 


« 




Discretized Integral Form 


The ultimate goal in computational fluids is to minimize the ap- 
proximations to the most fundamental modeling equations, presumably to 
salvage most of the physics, while attaining modest execution times on 
the available equipment. Hence the approach used here is founded on the 
discretization of the integral conservation law form of the Euler equa- 
tions, Eq. (3). This formulation, commonly referred to as a finite vol- 
ume method, yields the following discretized integral expression for a 
three-dimensional computational space with finite volume (cell) centers 
denoted i,j,k: 


where 
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In this expression the components of the dependent variable vec- 
tors, Q.- i ■_ represent average cell values. It is therefore evident 

1 > J 9 K 

that some method must be devised to accurately represent the vector-val- 
ued flux functions F, G, and H on the bounding surfaces (faces) of the 
cell. One of the methods used in this study is based on a one-dimen- 
sional analysis of the Riemann problem local to each interface, estab- 
lished by the discontinuous nature of the dependent variable vector Q in 
a finite volume field. To facilitate the introduction of this method, a 
return to a one-dimensional Cartesian space is in order. 
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First Order Flux Formula 


The discretized integral form of the Euler equations written for 
one spatial dimension appears as 
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Godunov^ proposed a procedure to obtain a global solution to Eq. ( 9 ) by 
solving the set of Riemann problems presented by the interface disconti- 
nuities. The Riemann problem local to each interface is costly to solve 
exactly, due to the necessary iteration. Many investigators'^ ’ 1 ® 
have made attempts to lessen this computational expense by approximating 
the solution of the Riemann problem. In essence these methods yield an 
approximate solution to the exact equation, and are referred to as ap- 
proximate Riemann solvers. 

In [11], Philip Roe suggests an alternate procedural choice. Roe 
proposed to obtain the exact solution to an approximate equation. The 
cleverness of Roe is evidenced by his choice of approximate equation. 
Consider the quasilinear form of Eq. (9)'s parent conservation law. 



where A(q^,qp) = 3f/3q is a constant matrix representative of local in- 
terface conditions. Matrix A is chosen to have the following specific 
list of properties, which Roe "christened Property U (since it is in- 
tended to ensure uniform validity across discontinuities)": 
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(i) It constitutes a linear mapping from the vector space q to the 
vector space f. 

— — — jjf 

(ii) As q L -*• q R -*• q, A(q L ,q R ) -*• A(q), where A = 

(iii) For any q L , q R , A(q L ,q R ) • (q R - q L ) = f R - f L . 

(iv) The eigenvectors of A are linearly independent. 

Restricting A to the satisfaction of Property U results in a spe- 
cial averaging process for the dependent variables from which A is con- 
structed. Referred to as "Roe averaged", the dependent variables are 
given by the following expressions: 
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where the total enthalpy, H, is defined 

H = - (e + p). (lid) 

P 

The flux difference can be expressed as 


df = f R - f L = A • (q R - q L ) = A • dq (12) 

where A is constructed with "Roe averaged" variables. Armed with the 
eigensystem of A and the knowledge that the interface differential dq is 
a proportional to the right eigenvectors of A, the interface flux dif- 
ference can be written relative to the right eigenvector basis as 


11 



df = l 


a.A (j) r (J) 

J 


= l 


+ a.A (j) r (j) 
J 


V» <J) r (J) 

J 


= df + df 


(13) 


Physically, the flux difference is shown to be the composition of a col- 
lection of waves. In this equation r^ is a right eigenvector of A; otj 
is the strength of the jth wave (the jump in the characteristic variable 
across it); A^ is an eigenvalue of A (the speed of the jth wave), and 
I - and E + denote summation over the negative and positive wave speeds, 
respectively. 

The interface flux, see Fig. (1), can be computed from either of 
the following formula: 
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This first-order interface flux formula developed for the one dimen- 
sional equations can be used in a multidimensional space provided the 
assumption that all waves travel normal to the interfaces is made. 


Higher Order (TVD) Flux Formulas 


In order to provide solutions of higher spatial accuracy, a family 

1 p 

of schemes can be represented by the addition of a corrective flux to 
the first order flux, Eq.(l4), produced by the preceding analysis: 
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The principle part of the truncation error for this flux formula is 
found to be 
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The details (common names, order of accuracy, etc.) of the members of 
this family can be found in [13,14], Two members, third-order (4> = 1/3) 
and fully upwind second order (<p = -1 ) were, by choice, exclusively used 
in this project. 

The. discussion of a higher order scheme inherently involves a 
method used to control spurious oscillations, i.e. dispersive errors. 
The method, actually methods, used in this project concern limiting com- 
ponents of the interface flux to produce total variation diminishing 
(TVD) schemes, i.e. non-oscillatory schemes. The following formulas 
have theoretical development as TVD schemes only in scalar nonlinear 
equations and systems of linear equations in one-dimension. One of the 
"limiters" used here, referred to as a minmod limiter, yields the fol- 
lowing expressions for the corrective flux terms: 
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where 


Ljtt.n) = minmod(a[^ /2 , Bo^ /2 ), (18) 

with a j ^ , a parameter proportional to the change in characteristic 
variables across nearby interfaces, defined by 
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where is a left eigenvector of A. The minmod limiter is then de- 

fined 


minmod[x,y] = sign(x) max{0, min[|x|, y sign(x)]} (20a) 


and the parameter B is a "compression" parameter given by 

1 < B < (20b) 

' -<P 

For this project the maximum B was used in all cases. 

Another limiter used here, this one credited to Roe 1 ^, is called 
Superbee. It is defined 
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where 

cmplim[x,y] = sign(x) max{0, min[jx|, By sign(x)], min[g|x| f y sign(x)]} 

(22) 

and the compression parameter B, not defined by Eq.(l8b), is taken here 
to be 2. 

Approximately Factored Implicit Scheme 


Up to this point no mention has been made as to what time level the 
numerical interface fluxes appearing in Eq. (8) are evaluated. The un- 
derlying theory of the approximate Riemann solver presented thus far is 
based on explicit concepts which result in an unattractive, rather 
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stringent time-step restriction. Equation (8) can be written in a 
linearized discrete- integral delta form to cover a broad class of im- 
plicit schemes 1 
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Some of the time differencing schemes represented are (9=1, \J>=1/2) three 
point backward, (6=1, 0 ) backward Euler, and (0=1/2, \J>=0) trapezoidal. 

Formally, all terms appearing in this equation should result from. a 
single theory. Superior results have been obtained, though, by evaluat- 
ing the residual term R n with flux difference split theory, and the 
left-hand-side (LHS) operator with flux vector split theory, see [17]. 
The rationale behind this is unclear and needs further investigation. 
Hence the following expressions complete Eq.(23) for this hybrid scheme. 
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( 25 ) 


R n = 6 .F + 6 .G + 6 .H, 
l J k 


where, F, G, and H result from flux difference split theory. 

The LHS of Eq. (23) tends to be cumbersome and difficult to invert, 
not to mention very costly. In light of this, the LHS was approximately 
factored into the product of two operators, each of which involve the 
passage of selected information. Here a forward and a backward operator 
are used as in [1] (block LU factorization), yielding the following two 
step (LU) scheme: 


[I 


8At 
1 +ip 


M- ] [ I 


6At 

1+ip 


M~ ]AQ 


At 

1+t|> 


l+i p 


AQ 


n -1 


(26) 


or 


[i + (SjA- + 5 B- + 6 r C-)]aQ 


Al_ R n + Jl AQ n -i 

1 +4> 1+t 


( 27 a) 


[i + ( 6 .A- + S.B~ + <5 k C-)]AQ n = AQ 


(27b) 


Q n+1 = Q n .+ AQ n (270 

Although factoring has been shown to degrade the unconditional sta- 
bility of Eq. (23), the (LU) scheme apparently retains this touted at- 
tribute. Equations (27) are in the final form of the mathematical model 
developed for the time-accurate analysis of turbomachinery. 

Since the approximate Riemann solver is a characteristic based ' 
scheme, the character istic variable boundary conditions developed in [ 3 ] 
relative to a three-dimensional time-dependent body-fitted reference 
frame for inflow, outflow, and impermeable boundaries are employed where 
applicable. As in [ 3 ] phantom cells are utilized to implement these 
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boundary conditions. The change in dependent variables, AQ n , is set to 
zero in the phantom cells for inflow, outflow, and impermeable bounda- 
ries . 
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LOGICAL APPROACH WITH COMPUTATIONAL OPTIMIZATION 
Selected Similarity Multi-Block 

As mentioned in the introduction, the scope of the field of turbo- 
machinery can encompass extremely complex flowfields. In order to ade- 
quately resolve these, a vast number of computational cells must be 
used. Several reasons can be cited for the segmentation of one virtu- 
ally insurmountable flow environment into several smaller more manage- 
able flow environments. This is commonly referred to as compositely 
gridding the.. field. There exist three basic methods of compositely 
gridding a field: overlaid, patched, and blocked. The approach taken 
here, a blocked grid method, is similar to that taken by Belk\ Belk 
has investigated many of the pitfalls encountered when constructing a 
code of a general multi-block nature. The emphasis in [M] was on devel- 
oping a computer algorithm to handle a completely arbitrary block ar- 
rangement, unfortunately this added to the complexity of the code. 

For the case of turbomachinery, the nature of the geometry suggests 
possible block arrangement and characteristic restrictions which can 
yield significantly simpler algorithm logic. The block structure pro- 
posed here for the specific case of dynamic cylindrical geometries (gen- 
erally, bladed or finned bodies of revolution) is referred to as 
selected similarity multi-block . 

Selected similarity multi-block is a method by which the benefits 
of a composite environment are attained while maintaining a reasonably 
simple, general algorithm. This is accomplished by requiring a fixed- 
map block arrangement with selectively-similar block characteristics. 
Hence by establishing prearranged, selectively-similar blocks with a 
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predefined solution path, algorithm logic can be simplified, thus en- 
hancing code readability. In short, an attempt is being made to walk 
the fine line separating code generality from code simplicity. The fun- 
damental restrictions applied in this approach deal with block arrange- 
ment and characteristics. 

A natural arrangement presents itself when dealing with cylindrical 
frames. Consider partitioning the macro-block (with accompanying global 
curvilinear axes) into partial cylinder blocks, see Fig. 2. Further 
considerations into the particular eccentricities of rotating geometries 
indicate a desire for additional partitioning in the axial direction, as 
shown in Fig. 3. Each partition in the axial, £, direction will hence- 
forth be referred to as a blade row, although it is not necessary for 
each axial partition to contain blades. Each blade row is granted lim- 
ited rotational freedom about the £ axis relative to the adjacent blade 
rows. With partitioning completed, a block referencing scheme must be 
installed. An obvious and simple choice is one which keys off the 
(global) curvilinear axes, also shown in Fig. 3. 

With an arrangement selected, the second fundamental restriction, 
selectively-similar block characteristics, can be addressed. To begin, 
the local block axes are assumed to follow the global axes relative to 
orientation and direction, i.e. there exists global similarity (see Fig. 
(M)). The block £,n,S (i,j,k) index limits (NI,NJ,NK) are restricted 
such that for all blocks within a blade row NI , NJ, and NK remain con- 
stant. In addition, NJ must remain constant between blade rows. NI and 
NK are allowed to vary between blade rows, although the blade row cir- 
cumferential cell count (BRCCC) must match. 
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For example consider Fig. (5), a three blade row configuration con- 
taining four partitions in the first blade row, i.e. four blade row 
blocks. Each of these inter-blade row blocks contain NI 1 xNJ 1 xNK 1 cells, 
where the subscript indicates the blade row. From the aforementioned 
restrictions on NJ, it is immediately established that NJ^N^^NJi . The 
BRCCC for blade row 1 is the product of NK 1 and the number of blade row 
1 blocks. Hence the BRCCC matching restriction implies that NK 2 and the 
number of blade row 2 blocks are limited to being factors of the BRCCC. 
In addition their product must equal the BRCCC. A similar condition 
exists for NK^ and the number of blade row 3 blocks. Summarizing these 
restrictions in tabular form 


Block Characteristic Restrictions 


INDEX 

SIMILARITY 

LIMITATION 

or I 

inter-blade row 

none 

or J 

global 

none 

or K 

inter-blade row 

BRCCC matching 


Memory Management 

Those. researchers with access to the genuine state-of-the art Cray 
2 with- its 256M words of internal memory might not concern themselves 
with in-core requirements. The vast majority of research though, is 
performed on equipment with lesser in-core capabilities. Also, it is 
understood that even Cray 2 users could find geometries that would push 
in-core requirements to the limits and beyond. Whence came an impetus 
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to manage the use of secondary memory, whether it is Cray’s rapid access 
Solid State Storage Device or any sequential storage medium for that 
matter. 

The use of computational blocks of differing cell count, and conse- 
quently differing memory requirement, necessitates the development of 
some special techniques to optimize memory utilization. A method re- 
ferred to as ribbon vector dynamic memory management^ is used here. As 
pointed out in [4], the principle of ribbon vector storage is frequently 
used to compress in-core storage to avoid wasting memory cells. The idea 
is to store the entire field on secondary memory with the exception of 
the block currently under execution. Ribbon vector dynamic memory man- 
agement allows the adjustment of the in-core field length to accommodate 
each block, regardless of size (limited, of course, by the total in-core 
capability), without excess memory cells. A version similar to that 
proposed by Belk has been integrated into this computer algorithm. 

Related to the topic of efficient memory use is the optimization of 
the processing of the memory. The time spent manipulating memory in- 
volves two general classes referred to as 10 wait time and CPU execution 
time. An effort to minimize both involves maximizing the efficient use 
of the hardware capability of the available equipment, as the ultimate 
limiting factor is the equipment itself. 

In the interest of minimizing 10 wait time, the use of unblocked 
data transmission to and from secondary memory, namely SSD, made signif- 
icant reductions. Essentially, all data contained within the ribbon 
vector is transferred between in-core and secondary memory via unblocked 
standard FORTRAN 10 statements. Some padding is required to maintain 
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proper unblocked data length (multiples of 512 words). Regarding CPU 
execution time, efforts are made to utilize vector hardware to its full- 
est, including a procedure to vectorize a backward or forward substitu- 
tion. 

The (LU) scheme involves point simultaneous solutions and backward 

(or forward) substitutions. The substitution portion poses quite a 

problem when attempting to utilize vector hardware. The difficulty is 

1 8 

referred to as a vector dependency or recursion problem. A method has 
been developed in order to circumvent this problem in a multi-dimen- 

p 

sional space. The procedure in a three-dimensional space is to "si- 
multaneously" process those cells lying on a special diagonal plane in 
computational space. Cells whose indices satisfy the equation of this 
diagonal plane, 

i + j + k = k (28) 

where k is a constant designating the plane level, are computationally 
independent. Hence with the aid of indirect addressing, diagonal plane 
processing is used to facilitate the vectorization of a backward or for- 
ward substitution. 

Interface Control 

Advancing a flowfield simulation in time requires solving the sys- 
tem model, Eqs. (27), at each computational cell within the global do- 
main for the new time level dependent variable vector Q. The solution 
path advocated here requires two sweeps through the global domain, one 
forward and one backward. The global path indicating the order in which 
the blocks are swept to insure analogous single block solution advance- 
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ment is shown in Fig. (6). All blocks are forward swept then all are 
backward swept. The forward sweep within each block consists of operat- 
ing on (applying Eq. (27a) to) each cell on a diagonal plane, beginning 
with the lowest level interior plane and advancing (forward) to higher 
level planes. The forward sweep is complete when the highest level in- 
terior plane has been updated. The backward sweep is similar although 
the planes are traversed in decreasing (backward) order, and Eqs. 
(27b, c) is the applied operator. 

Basic to the concept of multi-block is the requirement of continu- 
ous grid lines between blocks. This provides a relatively simple means 
of communication between blocks with stationary interfaces in computa- 
tional space. Values from within the domain of one block can be ex- 
tracted and then injected as phantom values in an adjacent block. As 
previously mentioned, each blade row can rotate relative to the adjacent 
blade rows. This rotation is presently governed by requiring continuous 
grid lines across the shearing block interface at all time levels. This 
is accomplished by maintaining equidistant circumferential cell spacing 
on the blade row boundaries. A shearing block interface in physical 
space yields a dynamic interface in computational space, i.e., the 
blocks continue to change communication partners. This complicates the 
simple communication procedure prescribed by a multi-block technique, 
but with appropriate index bookkeeping the extraction-injection proce- 
dure is still valid. 

When dealing with multi-block techniques, an issue of the accuracy 
necessary for the communication between blocks arises. Ultimately, the 
best communication possible is that analogous to a singly gridded field. 
The approach taken here with regard to internal block boundary communi- 
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cation is to mimic the cell communication within the blocks, thereby 
incurring no boundary induced error. Hence internal block boundaries 
are maintained with up to third-order spatial accuracy. For the case of 
re-entrant block boundary communication, periodicity and previous AQ or 
zero AQ are satisfactory. Block boundaries involving the farfield or 
impermeable surfaces are maintained as mentioned in the third section 
dealing with boundary conditions. 
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RESULTS 


Even Blade Count (GE UDF 8-8) 

The computer algorithm resulting from the procedures outlined in 
the previous sections has application in the simulation of even or une- 
ven blade count, single rotating, counter-rotating, or rotor-stator , 
axi-symmetric or non-symmetric, multi-stage, at angle of attack geome- 
tries. To support this claim several configurations are presently being 
analyzed, the results of which shall appear in subsequent publications. 
A configuration previously studied using a similar algorithm 1 ^ will be 
examined here. Since in [19] the interface fluxes were given by flux 
vector split theory and the block interfaces were only maintained to 
first order, it will be of interest to compare the two algorithms. 

The configuration is that of a counter-rotating unducted fan (UDF) 
emersed in an oncoming M=.72 axial flow, see Fig. 7. The configuration 
has two fan rows with eight blades per row. The highly swept, tapered, 
twisted, thin blades are designed to reduce the axial Mach number 
through the blading to alleviate compressibility losses. The forward 
row rotates clockwise and the aft row counter-clockwise. Both blades 
rotate with an advance ratio of 2.8. 

The solution appearing herein was obtained using only two blocks 
(benefitting from solution symmetry). Although only two blocks were 
used, axial interblock communication was implemented with a full buffer 
disk (temporary storage area for injected or extracted data) . The pro- 
cedure involves extracting data, imaging the data to form a full 360 
degree communication buffer disk, then allowing the appropriate data to 
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be injected based on the positional relationship between the blocks and 
the buffer disk. The mesh was H-type in all directions and each block 
contained 56x21x10 (i,j,k) cells. 

The solution (scalar quantities) of the UDF are periodic with re- 
spect to blade passages, i.e. the solution repeats each time the blades 
are re-aligned. Although this is the case, the velocity vector compo- 
nents in a Cartesian reference frame are only periodic with respect to 
blade position. The periodic behavior of scalar quantities gives a good 
measure of the "convergence" of such geometries. The convergence his- 
tory of the change in density, chosen to be defined as 
Log 1 0 ( (Ap )p ms / (Ap )p ms ] , is shown in Fig. 8. The transients associated 
with the arbitrary initial condition of uniform flow have diminished 
considerably after 3 complete revolutions of each blade row (480 time- 
steps), yet the results presented herein still retain some of the more 
stubborn transients. 

There exists evidence of the inherent unsteady behavior of the 
flow, though it is not prominent. Figures 9 and 10 support the comments 
in [19] regarding the lack of variation in the blade surface relative 
Mach number for both forward and aft blade rows. Though little varia- 
tion is shown overall, there is more variation in the aft blade row, as 
expected. Note, the previous comments regarding equidistant circumfer- 
ential cell spacing result in ten timesteps per blade passage. Hence 
each blade passage can be marked accordingly, with Step 1 representing 
an axial alignment of the blades (or as close to an alignment as one can 
get with highly swept, tapered, and twisted blades). Thus the suction 
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surface and pressure surface variations in Figs. 9 and 10 represent the 
respective maxima, and the steps indicated .are those from which these 
maxima were computed. 

Typical performance parameters are presented for comparison in 
Figs. 11. At first glance, one can see in Fig. 11a that the power coef- 
ficient (Power In) for the FDS scheme falls closer to experiment. The 
question arises whether this is the expected behavior. Consider the 
following plausible explanation. The FDS scheme is touted for less dis- 
sipation relative to the FVS scheme. Hence a reduction in numerical 
dissipation would tend to make for a 'more inviscid' flow, resulting in 
stronger, crisper shocks. This in-turn could yield pressure distribu- 
tions which produce lower torque coefficient values, hence lowering the 
power coefficient from that of the FVS scheme. Figures 12 and 13, indi- 
cating the blade surface relative Mach numbers for the two numerical 
algorithms, do not necessarily support this claim. It is believed the 
evidence needed to support this claim lies between the root and the mid- 
span plots, although since these did not appear in [19] a comparison 
cannot be made here. 

The reduced efficiency relative to the FVS scheme plotted in Fig. 
11b could also come from altered pressure distributions due to stronger 
shockwaves. In addition, the ratio of the blade torques given in Fig. 
11c indicate an increase relative to the FVS. It is difficult to access 
the cause of this increase, due in part to the influence on pertinent 
parameters such as flow velocity and flow angle relative to aft blade 
produced by the upstream blade row. Without benefit of a detailed flow- 
field snapshot between blade rows for both the FVS and the FDS schemes, 
no plausible arguments can be proposed. The conclusions drawn here re- 
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suit from a limited source of comparable data, further investigation is 
warranted. 
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CONCLUSIONS 


A computer program has been developed whose framework centers 
around the analysis of general rotating machinery. This general appli- 
cability has not been at the expense of algorithm simplicity. Painstak- 
ing thought has yielded an algorithm capable of simulating even or 
uneven blade count, single rotating, counter-rotat ing or rotor-stator , 
axi-symmetr ic or non-symmetric , multi-stage, at angle-of-attack geom- 
etries. The code presently stands at under 2000 lines, most of which 
are the numerics. A version of the numerical routine is resident on the 
NASA Marshall EADS Computer System. 

Future emphasis will focus on the development of simple axial in- 
terface control allowing non-uniform circumferential cell spacing. This 
could relax (or remove) the inherent dependency of the time-step on the 
dynamic interfaces. Parametric studies of the effect of time-step on 
solution accuracy, evaluating the compatibility of character istic vari- 
able boundary conditions with an approximate Riemann solver, and devel- 
opment of better grids for the simulated flowfields are in order. 
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Figure 5. Typical Configuration 




















Figure ,7. U DF Configuration 
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Figure 9. Forward Blade Row 
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Figure 10. Aft Blade Row 
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Figure 10. Aft Blade Row 
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Figure 11. Performance Parameters 
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Figure 13. Aft Blade Row 
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